################ Packages ##############################################
library(tidyverse)
library(fixest)
library(xtable)
library(ggpubr)
library(countrycode)
### functions
# IHS Transformation
ihs=function(x) log(x+sqrt(1+x^2)) # to revers ihs use sinh
# Demeans and IHS Transforms
idm=function(x) ihs(x)-mean(ihs(x),na.rm=T)
# Demeans
dm=function(x) x-mean(x,na.rm=T)
# Inverse Simpson
ihhi = function(x) 1/sum((x/sum(x))^2) 

################################### Set Working Directory ##########################################

# Getting the path of your current file
setwd(dirname(rstudioapi::getSourceEditorContext()$path))


################################### Data ##########################################
### main data
up = read.csv("../Data/data_main_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(stock_id)%>% # for demeaning at the stock level
  arrange(year)%>% # for calculating the lagged variables
  mutate(
    ffmsy = ihs(ffmsy), # Use only IHS because of the simulation (below). The demeaning is only important for the interacted variables, for the non-iteracted variables nothing changes
    gdp=idm(gdp_const),
    gdp2=idm(gdp_const)^2,
    tfp=idm(rtfpna),
    fao_ag_credit = idm(fao_agriculture_credit),
    fao_to_credit = idm(fao_total_credit),
    credit=idm(domestic_credit),
    credit_instrument = idm(credit_instrument),
    interest_factor=dm(1/(1 + lending_interest/100)),
    interest_factor = ifelse(inflation_cp<=30,interest_factor,NA), # This removes observations during hyper inflation
    inflation_factor=dm(1/(1+inflation_cp/100)),
    inflation_factor=ifelse(inflation_cp<=30,inflation_factor,NA), # This removes observations during hyper inflation
    bbmsy1 = idm(lag(bbmsy,1))
  )%>%
  as.data.frame()


### ram data
ram = read.csv("../Data/data_ram_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(country, stock_id, scientific_name)%>%
  arrange(year)%>% # for calculating the lagged variables
  mutate(ffmsy = ihs(ffmsy), # Fishing mortality
         gdp=idm(gdp_const), # GDP (World Bank)
         gdp2=idm(gdp_const)^2, 
         tfp = idm(rtfpna), # Total factor productivity (Penn)
         fao_ag_credit = idm(fao_agriculture_credit), # Credit to agriculture, fisheries, forestry from the FAO
         fao_to_credit = idm(fao_total_credit), # Total credit from the FAO
         credit=idm(domestic_credit), # Domestic credit to private sector (% of GDP)
         interest=idm(lending_interest), # Lending interest IMF
         inflation=idm(inflation_cp), # Inflation from World Bank
         interest_factor=dm(1/(1 + lending_interest/100)),
         inflation_factor=dm(1/(1+inflation_cp/100)),
         bbmsy1 = idm(lag(bbmsy,1)), # Fish stocks 
         species_cat_id =as.numeric(as.factor(isscaap_group)) # Species group
  )%>%
  as.data.frame()


# credit instrument data
za = read.csv("../Data/iv-shocks_final.csv")%>%
  mutate(residuals=feols(ihs(financial_export)~1|year + counterparty_country,data = .)$residuals)

zb = read.csv("../Data/iv-weights_final.csv")

################ Replicates Figure A.7: Shock distribution ##################################################
s1 = za%>%
  group_by(reporting_country,year)%>%
  summarise(financial_export = mean(financial_export,na.rm=T))%>%
  ungroup()%>%
  ggplot(aes(x = ihs(financial_export)))+
  geom_histogram(bins = 100)+
  theme_bw()+
  theme(text = element_text(family = "serif",size = 12))+
  labs(x = "Shock", y = "Count")

s2 = za%>%
  ggplot(aes(x = ihs(financial_export)))+
  geom_histogram(bins = 100)+
  theme_bw()+
  theme(text = element_text(family = "serif",size = 12))+
  labs(x = "Shock", y = "Count")

s3 = za%>%
  ggplot(aes(x = residuals))+
  geom_histogram(bins = 100)+
  theme_bw()+
  theme(text = element_text(family = "serif",size = 12))+
  labs(x = "Shock", y = "Count")



s = ggarrange(s1,s2,s3, labels = c("A)", "B)", "C)"),nrow = 1,
              font.label = list(face =  "plain", family ="serif",size=14)); s

ggsave(s,filename = "../Figures/iv_shock_distribution.pdf", device = "pdf",width = 16*1.5, height = 9, units ="cm")

################ Reproduces Figure A.8: Distribution of the A) average and B) bilateral weights ##################################################
w1 = zb%>%
  ggplot(aes(x = weight))+
  geom_histogram(bins = 100)+
  theme_bw()+
  theme(text = element_text(family = "serif",size = 12))+
  labs(x = "Weight", y = "Count")

zba = zb%>% # this intermediate step is needed for the summary statistics (see below)
  group_by(reporting_country)%>%
  summarise(weight = mean(weight,na.rm=T))

w2 = zba%>%ggplot(aes(x = weight))+
  geom_histogram()+
  theme_bw()+
  theme(text = element_text(family = "serif",size = 12))+
  labs(x = "Weight", y = "Count"); w2



w = ggarrange(w2,w1, labels = c("A)", "B)", "C)"),nrow = 1,
              font.label = list(face =  "plain", family ="serif",size=14)); w

ggsave(w,filename = "../Figures/iv_weight_distribution.pdf", device = "pdf",width = 16*1.5, height = 9, units ="cm")




################ Regression for Table A.2: Falsification test ##################################################

### the following code calculates base level and trend gdp and ffmsy
# base ffmsy in 1976 (the credit instrument data start in 1977)
ua = up%>%
  filter(year == 1976)%>% # 1977 is the first year of the credit data
  group_by(country)%>%
  summarize(base_ffmsy = mean(ffmsy,na.rm=T))

# pre ffmsy 10 year before the credit instrument
ub = up%>%
  filter(year == 1967)%>% # 1977 is the first year of the credit data
  group_by(country)%>%
  summarize(pre_ffmsy = mean(ffmsy,na.rm=T))

# same for gdp
ea = up%>%
  filter(year == 1976)%>% # 1977 is the first year of the credit data
  group_by(country)%>%
  summarize(base_gdp = mean(gdp_const,na.rm=T))

# same for gdp
eb = up%>%
  filter(year == 1967)%>% # 1977 is the first year of the credit data
  group_by(country)%>%
  summarize(pre_gdp = mean(gdp_const,na.rm=T))


## combines data

u = full_join(up,ua)%>%
  full_join(.,ub)%>%
  full_join(.,ea)%>%
  full_join(.,eb)%>%
  # calculates the trend between pre and base
  mutate(trend_ffmsy = (ihs(base_ffmsy) - ihs(pre_ffmsy)),
         trend_gdp = (ihs(base_gdp) - ihs(pre_gdp)))%>%
  filter(!is.na(base_ffmsy) & !is.na(base_gdp))%>%
  select(country,year,base_ffmsy,trend_ffmsy,base_gdp,trend_gdp,credit_instrument)%>%
  distinct()






r1 = feols(ihs(base_ffmsy) ~ ihs(credit_instrument)
           |year,
           cluster = ~ country + year, 
           data = u)

summary(r1)



r2 = feols(trend_ffmsy ~ ihs(credit_instrument)
           |year,
           cluster = ~ country + year, 
           data = u)

summary(r2)


r3 = feols(ihs(base_gdp) ~ ihs(credit_instrument)
           |year,
           cluster = ~ country + year, 
           data = u)

summary(r3)

r4 = feols(ihs(trend_gdp) ~ ihs(credit_instrument)
           |year,
           cluster = ~ country + year, 
           data = u)

summary(r4)




etable(r1,r2,r3,r4,
       digits = 3,
       digits.stats = 3,
       replace =T,
       dict=c("ihs(credit_instrument)" ="Instrument",
              "year" = "Year"),
       headers = c("Baseline extraction","Baseline extraction trend", "Baseline GDP","Baseline GDP trend"),
       depvar = F,
       file = "../Results/falsification_test.tex")


############################ Replication for Table A.3: IV excluding countries with high credit dependency on individual countries ##############
### this code excludes all countries with at least one bilateral weight greater than 0.5
hwc = zb%>%
  group_by(country = counterparty_country)%>%
  summarise(max_weight = max(weight))%>%
  filter(max_weight<0.5)%>%
  as.data.frame()
  

# replicates the main specification from the main results
h1 = feols(ffmsy ~  pr 
           + gdp + gdp2 
           | stock_id + country + year
           |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
           cluster = ~ country + year + scientific_name, 
           data = up)


summary(h1)

# replicates the main specification but without high-weight countries
h2 = feols(ffmsy ~  pr 
           + gdp + gdp2 
           | stock_id + country + year
           |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%
             filter(country%in%hwc$country))


summary(h2)

etable(h1,h2,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit"),
       dict=c("pr" ="Property rights",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "credit x pr",
              "credit:pr" = "Credit x Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       headers = c("Main","Excluding coutries with large weights"),
       depvar = F,
       file = "../Results/high_weight.tex")



